bellman.m
Calculates LIFO equilibrium continuation values using Bellman iteration and plots the results.
Contents
Set up the problem's parameters.
This section first checks to see if the parameters of the problem are defined in the structure bellmanARG . If this does not exist, then it defines them.
if exist('bellmanARG','var') %Parameters for constructing the Markov chain. markovARG=bellmanARG.markovARG; %Producer's surplus per customer as a function of N. piF=bellmanARG.piF; %Sunk entry cost and per-period fixed cost. phi=bellmanARG.phi; kappa=bellmanARG.kappa; %Discount factor beta=bellmanARG.beta; else %Parameters for approximating the innovation distribution approximateARG.F=@(x) (1+erf(x/sqrt(2)))./2; approximateARG.Frange=3.0; approximateARG.k=51; %Parameters for constructing the Markov chain. markovARG.approximateARG=approximateARG; markovARG.rho=1.0; markovARG.sigma=0.30; markovARG.omegaCenter=0.0; markovARG.omegaStep=0.01; markovARG.omegaWidth=3; %Parameters describing profits and %discounting. k=4; piF = @(x) k*ones(size(x)); %Producer surplus per customer kappa=1.75; %per-period fixed cost beta=1.05^(-1); %5 percent annual interest rate. phi = @(N) 0.25*beta/(1-beta); %Sunk cost of entry end
Create Markov chain if required.
The results are stored in Pi and omega.
if ~exist('Pi','var') markov end
Calculate upper bound on firm numbers.
nCheck=find((max(exp(omega))*piF(1:1:100)./(1:1:100)>=kappa),1,'last');
Solve individual entrants' problems
- Use Bellman iteration
- Update the transition rule for N
% Store the individual firms' value functions in a three-dimensional array. % The final dimension gives the rank of the firm. % The second dimension gives the number of firms presently active. % The first dimension indexes the demand states. vFuncs=zeros(nCstates,nCheck,nCheck); for R=nCheck:-1:1 %Initialize the value function. v=zeros(nCstates,nCheck-R+1); %Initialize the transition rule for N if necessary. if R==nCheck; nPrime=nCheck*ones(nCstates,1); end %Create nPrimeIndex, which determines which column of v(:,:) contains the relevant continuation value. if R==nCheck nPrimeIndex=ones(nCstates,nCheck-R+1); else nPrimeIndex=zeros(nCstates,nCheck-R+1); for i=R:1:nCheck; nPrimeIndex=nPrimeIndex+(nPrime==i)*(1+(i-R)); end end % Bellman Iteration tol=1.0e-7; vDistance=1; while vDistance>=tol %Use nprime to assign the appropriate continuation value. Piv=Pi*v; Ev=zeros(nCstates,nCheck-R+1); for i=1:nCstates; Ev(i,:)=Piv(i,nPrimeIndex(i,:)); end vPrime=beta*(Pi*exp(omega')*ones(1,nCheck-R+1)).*piF(nPrime)./nPrime-beta*kappa+beta*Ev; %Take the maximum of the expected continuation value and zero. vPrime=max(vPrime,0); %Measure the distance between the initial and new value functions. vDistance=max(max(abs(vPrime-v))); %Update the value function. v=vPrime; end %Store the value function if R>1; %In this case, we need to pad the matrix of values with NaNs vv=[NaN(nCstates,R-1) v]; else vv=v; end vFuncs(:,:,R)=vv; %Update nprime for the problem of a firm with rank R-1 %Take away one firm if exit is optimal. nPrime=nPrime-(v==0); %Add a column to nPrime for the initial state N=R-1 (The Matlab editor complains about the speed consequences of %this assigment. For a big problem where the speed drag from copying nPrime to a new location in memory is too much %to pay, preallocate nPrime to a nCstates x nCheck+1 array and only use the relevant column in the Bellman equation %above.) nPrime=[(R-1)*ones(nCstates,1) nPrime]; %Add firms to the new column in the states where entry is optimal. for i=R:1:nCheck nPrime(:,1)=nPrime(:,1)+(squeeze(vFuncs(:,i,i))>phi(i)); end end
Create final version of nPrimeIndex.
nPrimeIndex=zeros(nCstates,nCheck+1); for i=1:1:nCheck; nPrimeIndex=nPrimeIndex+(nPrime==i)*(1+i); end
Determine the ergodic set for N.
minN=min(min(nPrime)); maxN=max((max(nPrime).*(max(nPrime)>(0:1:nCheck)))); %Highest value of N transited to from a lower value. %Trim the matrix nPrime to eliminate states outside of the ergodic set. nPrime=nPrime(:,minN+1:maxN+1); nPrimeIndex=nPrimeIndex(:,minN+1:maxN+1)-minN+1; %Trim the value functions to eliminate values %of N above the ergodic set. vFuncs=vFuncs(:,1:maxN,1:maxN);
Plot the value functions.
This section creates a raw inelegant plot which can be used to verify that the value functions satisfy the restrictions of theory.
if ~exist('bellmanARG','var') figure(1) set(gcf,'Units','inches','Position',[0 0 4 6]); for i=1:maxN; subplot(maxN,1,i); vi=squeeze(vFuncs(:,:,i)); plot(omega'*ones(1,maxN-i+1),vi(:,i:maxN)); %This line plots the value function's many ``branches'' all at one go. hold plot(omega',phi(i)*ones(size(omega))); %This line adds a horizontal line to mark the entry cost. end end
Current plot held Current plot held Current plot held Current plot held
